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^ Abstract 

i cqj ; 

Latin hypercube designs (LHDs) with space-filling properties are widely used for 
emulating computer simulators. Over the last three decades, a wide spectrum of LHDs 
have been proposed with space-filling criteria like minimum correlation among factors, 
maximin interpoint distance, and orthogonality among the factors via orthogonal arrays 
(OAs). Projective geometric structures like spreads, covers and stars of PG(p — l,q) 
can be used to characterize the randomization restriction of multistage factorial exper- 
iments. These geometric structures can also be used for constructing OAs and nearly 
^ ■ OAs (NOAs). In this paper, we present a new class of space-filling LHDs based on 

NOAs derived from stars of PG(p — 1,2). 
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1 Introduction 

Latin hypercube sampling scheme was first proposed by McKay et al. (1979) as an alternative 

to random sampling in the Monte Carlo methods for numerically integrating complex multi- 
• 

^ | dimensional functions. Later on, the Latin hypercube designs (LHDs) became very popular 
in computer experiments for building statistical metamodels (Santner et al. 2003). Ran- 
dom LHDs can easily be constructed and are also useful in some applications, for instance, 
numerical integrations. However, random LHDs are not always suitable from a modeling 
viewpoint, for example, if all points are aligned along the main diagonal of the input space. 

Space-filling designs have become very popular in computer experiments for building 
surrogates of the deterministic simulators. Since replicate runs of a deterministic simulator 
generate identical outputs, it is preferred that the inputs are spread out in the design space. 
Popular space-filling designs proposed thus far include, A, D, /-optimal designs (for minimiz- 
ing generalized variance), mean squared error based designs, distance based designs (maximin 
and minimax), uniform designs, space-filling LHDs, orthogonal arrays (OAs), scrambled nets 
and Sobol sequences (see Santner et al. (2003) and Fang et al. (2006) for an overview). 
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In this paper, we discuss space-filling LHDs, a popular class of designs in computer 
experiments. Over the last three decades, a wide spectrum of LHDs have been proposed with 
different space-filling criteria, for instance, minimum correlation among factors (Iman and 
Conover 1982), maximin interpoint distance (Morris and Mitchell 1995), and orthogonality 
among the factors via OAs (Owen 1992; Tang 1993). The construction of LHDs with space- 
filling criteria like maximin distance or minimum correlation often requires computationally 
intensive search, whereas, OA-based LHDs are easy to construct as long as the OAs exist. 

The existence of a desired OA is not always guaranteed, and the construction can also be 
challenging (Hedayat et al. 1999). OAs can be constructed using a variety of combinatorial 
objects like linear codes, difference schemes and mutually orthogonal Latin squares. Rains et 
al. (2000) discussed the existence and construction of OAs using a spread of a finite projective 
geometry, V = PG(p—l, q), and called them geometric OAs. A spread of V is a set of disjoint 
subspaces or flats of V that cover V . Ranjan et al. (2009) established the equivalence between 
2 P factorial experiments with multiple randomization restrictions and different geometric 
structures of PG(p—l, 2) (e.g., spreads and covers). Here, a randomization defining contrast 
subspace (RDCSS) with t(< p) independent randomization factors is equivalent to a (t— 1)- 
flat of PG(p — 1,2). This paper focuses on a class of LHDs that are based on arrays 
constructed using RDCSSs of a 2 P factorial experiment (i.e., we only consider q — 2). 

For efficient analysis of a multistage factorial experiment, it is desirable to construct 
disjoint RDCSSs. Ranjan et al. (2009) established the existence of a desired set of disjoint 
RDCSSs via the existence of a spread of a PG(p—l, 2). However, in many practical situations 
(e.g., the plutonium alloy experiment of Bingham et al. 2008), overlap among the RDCSSs 
cannot be avoided. For such cases, Ranjan et al. (2010) proposed designs based on a new 
geometric structure called a star - a set of distinct flats of PG(p — 1, q) that share a common 
overlap (the nucleus). That is, the proposed design is a set of distinct RDCSSs with a 
common overlap. A star that is also a cover (covering star) of PG(p — l,q) simplifies to a 
spread if the nucleus is empty. 

We have discovered a new class of space-filling LHDs that can be constructed using stars 
of PG(p — 1,2). It turns out that a star with non empty nucleus generates nearly OAs 
(NOAs). In the spirit of Rains et al. (2000), we refer to these star-based NOAs as geometric 
NOAs. By following Tang's OA-based LHD construction algorithm, we construct a class 
of geometric NOA-based LHDs. Although such LHDs are not always very space-filling, a 
near orthogonality (e.g., Xu and Wu 2001) or space-filling criterion can be used to search 
for a good one. To avoid the search, we propose a set of guidelines for carefully distributing 
the factorial effects among RDCSSs of the star which ensures space-filling LHDs. Moreover, 
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the existence of OA-based LHDs are limited to only few n x d combinations, whereas, the 
existence conditions for stars are less stringent. 

The remainder of the paper is organized as follows. Section 2 presents an overview of 
LHDs, RDCSSs in a TP factorial experiment, and spreads and stars of a PG(p— 1, q). In Sec- 
tion 3, we establish theoretical results for the existence and an algorithm of the construction 
of geometric-NOAs. Section 4 concludes the paper with a few remarks. 

2 Background 

This section starts with a brief review of random LHDs and OA-based LHDs. Then a few 
results are presented to establish the equivalence between a multistage factorial design with 
randomization restrictions and geometric structures of PG(p — 1,2). 

2.1 Latin hypercube designs 

Let L(n, d) be an LHD with n runs in d factors, where Lij denote the level of factor j in 
the i-th experimental run, and each factor includes n uniformly spaced levels. In computer 
experiments, the input spaces are typically bounded hyperrectangles, and can be transformed 
to unit hypercubes. A random L(n, d) in [0, l] d has = (vrj(i) — Uij)/n for 1 < j < d and 
1 < i < n, where Uij ~ Unif(0, 1) and {^(1), 7ij(n)} is a random permutation of {1, n} 
(see Tang (1993) for details). As a result, there are (n\) d distinct LHDs, even if we ignore 
the slight variation due to Unif(0, 1) perturbations. 

Although LHDs have nice one-dimensional projection property, that is, one point each 
in ((i — l)/n,i/n) for 1 < i < n, random LHDs can be quite undesirable from modeling 
viewpoint. Figure [1] presents two realizations of random LHDs in [0, l] 2 . 




0-2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 



(a) A good design for modeling (b) A bad design for modeling 

Figure 1: Two realizations of random LHDs in [0, l] 2 . 
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Space-filling LHDs are commonly used in computer experiments for building statistical 
surrogates of complex computer simulators. In this paper, we propose star-based LHDs 
which are generalizations of the OA-based LHDs. 

Definition 1 An N x d array A denoted by OA(N, S1S2 ■ ■ ■ Sd, r) is said to be a strength r 
OA with N runs and d factors, if factor j has Sj levels {0, Sj — 1} and each N x r subarray 
contains every possible r -tuple an equal number of times. 

The special case of si = s 2 = • • ■ = Sd corresponds to a symmetric OA denoted by 
OA(N, s,d,r). Although mixed-level (or asymmetric) OAs have been investigated in recent 
years (see Hedayat et al. (1999) for references), it is a relatively much less explored area as 
compared to the symmetric OAs. 

A simple existence condition of an asymmetric OA follows from the strength aspect of 
Definition [TJ that is, N must be a multiple of s^sJf ' ' ' s p d f° r every set of u\, Ud G {0, 1} 
such that Y2i=i u i — r - Another popular existence result comes from the Rao bound: N — l > 
Yli=i( s i ~ !)• Determining the existence of OAs is nontrivial, and becomes more and more 
difficult as the strength r and the number of levels Sj increase. Even OAs with strength 2 
do not always exist for arbitrary N and d (see Hedayat et al. (1999) and Rains et al. (2002) 
for more results). 

An OA(N, S1S2 ■ ■ ■ Sd, r)-based LHD will have iV runs and d factors with = (iTj(Aij) + 
Uij)/sj for 1 < j < p and 1 < i < N, where ity ~ Unif '(0,1) and iTj are independent 
random permutation of {0, ...,Sj — 1}. As a result, the input space [0, l] d is partitioned in 
to a si x S2 x • • • x Sd grid and N points are assigned to the cells according to the scheme 
directed by the rows of OA. Tang (1993) showed that LHDs based on a strength r OA 
inherit r-dimensional projection property. However, such LHDs do not guarantee projection 
property (i.e., space-filling) for lower dimensional hyperspace (see Figure [2]). 




0.0 0.2 0.4 0.6 0.8 1.0 



Figure 2: An OA(9, 4, 3, 2)-based LHD in [0, l] 2 . 
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In this paper, we propose a new class of star-based geometric NOAs for constructing 
space-filling LHDs. As mentioned earlier, OAs exist for only a few n x d combinations, 
whereas, the existence conditions for stars are less stringent. 

2.2 RDCSSs and Projective Geometries 

Ranjan et al. (2009) proposed a unified framework using finite projective geometry for 
the existence and construction of multistage factorial designs with randomization restriction 
(e.g., nested designs, split-plot designs, split-lot designs, and combinations thereof). Let 
V = PG(p — 1,2) be the (p — l)-dimensional finite projective space defined by the set of 
all p-dimensional pencils over GF(2). Here, a p-dimensional pencil with r nonzero elements 
uniquely corresponds to an r-factor interaction, and V denotes the set of all 2 P — 1 factorial 
effects (excluding the grand mean). 

The restrictions on the randomization of experimental runs lead to grouping experimental 
units into sets of trials. These sets are formed by using linearly independent pencils of 
V, also referred to as the randomization restriction factors (like the blocking factors in a 
blocked factorial design). A set S of all non-null pencils formed as linear combinations 
of t independent randomization restriction factors in V constitutes a (t — l)-dimensional 
subspace (or, (t — l)-flat) of V (here \S\ = 2* — 1). We call such a subspace a t-dimensional 
randomization defining contrast subspace (RDCSS). 

For efficient analysis of such experiments, it is desirable to construct disjoint RDCSSs. 
Ranjan et al. (2009) established the existence of a set of disjoint RDCSSs in a TP factorial 
experiment via the existence of a spread of a PG(p — 1,2). 

Definition 2 For 1 < t < p, a balanced (t — l)-spread of V = PG{p — l,q) is a set, ip, of 
(t — 1) -flats ofV which partitions V . 

The size of a balanced (t — l)-spread if> of PG(p — l,q) is \ip\ = (q p — l)/(<7* — 1). A 
necessary and sufficient condition for the existence of a (t — l)-spread is that t divides p 
(Andre 1954). See Ranjan et al. (2009) for more results on the balanced and mixed (partial) 
spreads. In many practical situations (e.g., the plutonium alloy experiment of Bingham et 
al. (2008)), the overlap among the RDCSSs cannot be avoided. For such cases, Ranjan et al. 
(2010) proposed designs based on a new geometric structure called a star - a set of distinct 
flats of PG(p — 1, q) that share a common overlap (the nucleus). 

Definition 3 A balanced covering star Q = St(/j,,t,to) of V = PG(p — l,q) is a set of \i 
rays ((t — 1)- flats) and a nucleus ((to — 1)- flat), where the nucleus is contained in each of 
the fx rays (i.e., t <t), and /i = (q p ~ to — 1) / — 1) . 
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A necessary and sufficient condition for the existence of a balanced covering stars St(p,, t, to) 
of PG(p—l, q) is (t — to) divides (p — to)- Let St(t 1 , t^; t ) be a mixed or unbalanced cover- 
ing star of PG(p — 1, q) with fi rays and (to — l)-dimensional nucleus, such that t\ < • • • < t^. 
The next two lemmas are taken from Ranjan et al. (2010). 

Lemma 1 For the existence of a covering star St(ti, t M ; to) of V — PG(p — l,q), the 
following conditions are necessary: 

(i) q P-to _ i = £f =1 (^-*> _ i), 
(ii) U + tj-to<p for every i ^ j = 1, k), 

A necessary and sufficient condition for unbalanced covering star is still unknown. The 
next result guarantees the existence of a balanced star for every t and p (t < p). 

Lemma 2 For every t (2 < t < p) and to — t — 1, there exists a balanced covering star 
St(ji,t,t ) ofV = PG(p-l,q) with fi = (qP- t+1 - l)/(q - 1). 

Next we discuss how these star-based multistage factorial designs can be used to construct 
space-filling LHDs. 

3 Star based LHD 

In this section, we first establish new existence results for geometric NOAs that are derived 
from stars. Then we propose an algorithm (a generalization of Rains et al. (2002)) for 
constructing such NOAs. Finally, we present a few guidelines for constructing NOAs that 
lead to space-filling LHDs. 

Theorem 1 The existence of a covering star St(ti, ...,t M ;to) of PG(p — 1, 2), is a sufficient 
condition for the existence of an NOA(2 p , s^s^ ■ ■ • Si k , 2) with l<k<fi,l<ii<i2< 
• • • < «fe < A* an d s j = 2 <J for all 1 < j < fi. 

Similarly, the existence of a balanced star St(fi;t;to) of PG(p — 1,2) suffices the ex- 
istence of NOA(2 p , 2', k, 2) for all 1 < k < \i. Moreover, it turns out that the maxi- 
mal NOA(2 p , 2*, fx, 2) is a 2~*° fraction of an OA(2 p+t °, 2*, /i, 2), and its existence requires 
(t — to) 10° — to) which follows from the existence of a balanced star. 

Theorem 2 For every 1 < t < p, there exist NOA(2 p , 2*, k, 2) for all I < k < 2 p ~ t+1 -l, and 
the maximal NO A{2 P ,2 l ,2 p - t+l -1,2) is a2~^ fraction of an OA(2 p+t -\ 2*, 2 p ~ t+1 - 1, 2). 
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Proofs of Theorems [T] and [2] follows from Lemma [T] and (H and the construction of such 
NOAs using covering stars of V — PG(p—l, 2). Algorithm [JJ is a generalization of the spread 
to OA construction recipe presented in Rains et al. (2002). Note that a covering star with 
empty nucleus (i.e., t = 0) is a spread. 

For this construction, we use the binary (vector or pencil) representation of all effects. 
Let {a ly a 2 , a 2 p-i} be the ordered set of all effects in V, and St(ti, ■■■,t fJi ; t ) be a covering 
star of V with /i rays {R\, ...,-R^}, where \Rj\ = 2 tj (= Sj, say). Algorithm Q] constructs an 
NOA(2 p , sis 2 • • • Sp, 2) denoted by A = [A* x : A* 2 : ■ • • : A^]. 

Algorithm 1 Star to NOA construction 
1: for j ; = 1 — > n do 

2: For j-th stage of randomization restriction, find tj independent randomization re- 
striction factors {5[ j \ ...,5$} such that R d = (5[ j) , ...,6$). 
3: for i = 1 -> 2 P - 1 do 
4: for Z = 1 — > tj do 

5: Compute b\f = a { ■ # inner product over mod(2) 

6: end for 

7: Define A^ = J2tAj2 tj ~ l - 

8: end for 
9: end for 



For every j G {1,2,..., /i}, the i-th element of the NOA is A {j e {0, 1, 2^ - 1} for all 
1 < i <2 P — 1. As a convention, we append a row of zeros at the beginning. We now present 
a few examples to illustrate the theoretical results and Algorithm [TJ 

Example 1 Let Q be a covering star of PG(3, 2) (or, in a 2 4 experiment) such that the 
rays (or RDCSSs) are of size seven each (i.e., p = 4 and t = 3). Here, {ai, ...,015} = 
{-D, C, CD, B, ABCD}. Since there always exists a covering star for t = t — 1 (with /x = 
2 4-2 — 1 = 3), Theorem [2] ensures the existence of NOA(16, 8, fc, 2) for all 1 < k < 3. Follow- 
ing the construction in Ranjan et al. (2010), let Q = {Ri, R2, R3} with Ri = (A,B,ACD), 
R 2 = (C, D, ABC) and R 3 = (AC, BC, AD). Note that the nucleus is (AB, CD). 

The first column of the maximal NOA is A*i = (0, An, A21, Ai5^) T , with, An = 
&S2 3 - 1 + 4?2 3 ~ 2 + 4? 23 " 3 - A s an example, for % = 2, a 2 = C and bf} = (1,0,0,0) T ■ 
(0, 0, 1, 0) = 0, 6g } = (0, 1, 0, 0) T • (0, 0, 1, 0) = and &g } = (1, 0, 1, if ■ (0, 0, 1, 0) = 1. 
Therefore, A21 = 0- 4 + 0- 2 + 1- 1 = 1. The transpose of the final NOA is as follows: 

/0110233254457667\ 
A T = 0257134613460257 
^0167234554327610 
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Though LHDs derived from star-based NOAs have nice geometric features, not all such 
LHDs are very space-filling. One can use a near orthogonality (Xu and Wu 2001) or space- 
filling criterion to choose a suitable NOA for constructing space-filling LHD, however, the 
search may still be computationally intensive. To avoid the search, we propose a few guide- 
lines (not exhaustive) for carefully distributing the factorial effects among the RDCSSs of 
the star which ensures space-filling LHDs. 

• (Gl) 81 (the first randomization factor in any ray) is not an element of the nucleus 

• (G2) 5 i l Jl) ^ 5\ n) for all I 

. (G3) 6<rt + 6W ? sW + Sjp for all (h,l 2 ) e {1,2, ...,**} x {1, 2, t h }} 

These guidelines suggest a careful selection of the randomization restriction factors for 
constructing the desirable geometric NOA. In the next few examples we use the star in 
Example 1 but choose ^ s carefully to illustrate the importance of these guidelines. 

Example 2 Consider the same star as in Example [H that is, a covering star Q = St(3; 3; 2) 
of PG(3,2) (with p = 4, t = 3 and t = 2), however, choose <f } 's such that only "Gl" 
is violated. Let R x = (AB, B, ACD), R 2 = (D,C,ABC) and R 3 = (AC,BC,CD). Note 
8i = AB belongs to the nucleus. Figure [3] shows the corresponding LHDs with identity 
permutation and no Unif(0, 1) perturbation. 



5 6 7 



5 6 7 



5 6 7 



(a) A bad design (i?i and R2) (b) A bad design and R3) (c) A good design (R 2 and R3) 
Figure 3: Three realizations of NOA(16, 8, 2, 2)-based LHDs in [0, l] 2 . 



It is clear from Figure |3] that LHDs derived using Ri (Figures 3(a) and 3(b)) are not as 



space-filling as compared to that obtained without it (Figure 3(c)). 



Example 3 Consider the same star as in Examples [j] and [2], however, choose <Jj s such 
that only "G2" is violated. Let i?i = (A,B,ABCD), R 2 = (C,D,ABCD) and R 3 = 



(AC,BD,BC). Note 5 : 



(i) 



6. 



(2) 



ABCD, which also belongs to the nucleus. Figure H] 



shows the corresponding LHDs with identity permutation and no Unif(0, 1) perturbation. 



5 6 7 



3 4 5 6 7 



5 6 7 



(a) A bad design (i?! and R 2 ) (b) A good design (Ri and R 3 ) (c) A good design (i? 2 and i? 3 ) 
Figure 4: Three realizations of JVOA(16, 8, 2, 2)-based LHDs in [0, l] 2 . 



It is clear from Figure 4(a) that LHD based on (Ri, R 2 ) shows correlation among the two 



factors (induced due to 6 



5. 



(2) 



ABCD), and may not be very desirable for modeling. 



LHDs shown in Figures 4(b) and 4(c) ) are reasonably space-filling in two-dimensional space. 



Example 4 Consider the same star as in Examples [T] - [3X however, choose 5\ ; s such 
that only "G3" is violated. Let R 1 = (A, B, ACD), R 2 = (C, ABD, ABC) and R 3 = 



(AC,AD,BC). Note that (Ri,R 3 ) violate "G3" as 8^ + 5 3 L> = B + ACD = ABCD and 



4 3) + 5 ( 3 ] = AD + BC = ABCD, and (R 2 , R 3 ) violate "G3" as 5\ z> + 5^ Zj = C + ABC = AB 



( 2 ) , X( 2 ) 



and 5^ + 5^ = AC + BC = AB. Figure [5] shows the corresponding LHDs. 



5 6 7 



3 4 5 6 7 



5 6 7 



(a) A good design (Ri and R2) (b) An okay design (Ri and R3) (c) A bad design (R2 and R3) 
Figure 5: Three realizations of NOA(l6, 8, 2, 2)-based LHDs in [0, l] 2 . 
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It is clear from Figure [5] that the violation of "G3" caused adverse effect on the random- 
ness of the points and space-filling behaviour of the LHD. 

Example 5 Consider the same star as in Examples [1] - [4j however, we wish to choose 
Si's that follow all three guidelines (i.e., no violations). Let R\ = (B,ACD,AB), R2 = 
(D,C,ABC) and R 3 = (AC,BC,CD). Note that the set {R u R 2 , R 3 } here is very similar 
to that in Example [21 except the order of 5 Z (randomization restriction factors for Ri). 
Figure [6] depicts the corresponding LHDs. 



5 6 7 



2 3 4 5 6 7 



5 6 7 



(a) A good design (Ri and R2) (b) A good design (Ri and R3) (c) A good design (R 2 and R3) 
Figure 6: Three realizations of NOA(l6, 8, 2, 2)-based LHDs in [0, l] 2 . 



It is clear from Figure that the NOA derived from this star generates LHD with two- 
dimensional space-filling projections. 



4 Concluding Remarks 

This paper presents a new class of space-filling LHDs that are based on geometric NOAs 
derived from covering stars - designs for multistage factorial experiments with random- 
ization restrictions. Though we assumed two-level factorial designs (i.e., PG(p — 1,2)) 
in Section 3, theoretical results and construction algorithms can easily be generalized for 
mixed/unbalanced covering or partial stars in PG(p — l,q). 

Although we have proposed a few guidelines for constructing space-filling NOA-based 
LHDs, it may not always be possible to satisfy all conditions (Gl, G2 and G3) in a star 
construction. In such a case, one can use a near orthogonality (Xu and Wu 2001) or space- 
filling criterion to choose a suitable NOA for constructing space- filling LHD. 

In Example [U the LHDs are based a covering star St(3; 3; 2) of PG(3, 2), i.e., we chose 
t = p — 1. One could instead use another 2 < t < p for constructing a non-trivial star. 
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For instance, t = 2 and to = t — 1 would generate a covering star St(7;2; 1) of PG(3,2). 
The advantage of using St(7; 2; 1) instead of 5^(3; 3; 2) would be the opportunity to construct 
LHDs with d = 7, however, the number of levels per factor will be reduced from eight to four. 
That is, the design would be slightly less space-filling in the one- dimensional hyperspace. 
Thus large covering stars are more useful. 
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